/*
 * NormalizedSteinharHartModelThermistor.cpp
 *
 * Copyright 2012 Pieter Agten
 *
 * This file is part of tltgen.
 *
 * tltgen is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * tltgen is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with tltgen.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * @file   NormalizedSteinhartHartModelThermistor.cpp
 * @author Pieter Agten (pieter.agten@gmail.com)
 * @date   2 sep 2012
 */


#include <boost/lexical_cast.hpp>
#include <cmath>
#include <stdexcept>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "ThermistorVisitor.hpp"

#include "NormalizedSteinhartHartModelThermistor.hpp"

namespace tltgen {

  NormalizedSteinhartHartModelThermistor::NormalizedSteinhartHartModelThermistor(Resistance r0, double a, double b, double c, double d) :
    r0(r0), a(a), b(b), c(c), d(d)
  {
    if (! isValidR0(r0)) {
      throw std::invalid_argument("Value '" + boost::lexical_cast<std::string>(r0) + "' is invalid for parameter R0.");
    }
    if (! isValidA(a)) {
      throw std::invalid_argument("Value '" + boost::lexical_cast<std::string>(a) + "' is invalid for parameter A.");
    }
    if (! isValidB(b)) {
      throw std::invalid_argument("Value '" + boost::lexical_cast<std::string>(b) + "' is invalid for parameter B.");
    }
    if (! isValidC(c)) {
      throw std::invalid_argument("Value '" + boost::lexical_cast<std::string>(c) + "' is invalid for parameter C.");
    }
    if (! isValidD(d)) {
      throw std::invalid_argument("Value '" + boost::lexical_cast<std::string>(d) + "' is invalid for parameter D.");
    }
  }

  NormalizedSteinhartHartModelThermistor::~NormalizedSteinhartHartModelThermistor() { };

  bool NormalizedSteinhartHartModelThermistor::isValidR0(Resistance r0)
  {
    return tltgen::isValidResistance(r0);
  }
  Resistance NormalizedSteinhartHartModelThermistor::getR0() const
  {
    return this->r0;
  }

  bool NormalizedSteinhartHartModelThermistor::isValidA(double a)
  {
    return std::isfinite(a) && a > 0;
  }
  double NormalizedSteinhartHartModelThermistor::getA(void) const
  {
    return this->a;
  }

  bool NormalizedSteinhartHartModelThermistor::isValidB(double b)
  {
    return std::isfinite(b) && b > 0;
  }
  double NormalizedSteinhartHartModelThermistor::getB(void) const
  {
    return this->b;
  }

  bool NormalizedSteinhartHartModelThermistor::isValidC(double c)
  {
    return std::isfinite(c);
  }
  double NormalizedSteinhartHartModelThermistor::getC(void) const
  {
    return this->c;
  }

  bool NormalizedSteinhartHartModelThermistor::isValidD(double d)
  {
    return std::isfinite(d);
  }
  double NormalizedSteinhartHartModelThermistor::getD(void) const
  {
    return this->d;
  }

  bool NormalizedSteinhartHartModelThermistor::isValidParameterIndex(unsigned int pi) const
  {
    return 0 <= pi && pi < 5;
  }

  double NormalizedSteinhartHartModelThermistor::getParameterValue(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument("The given parameter index is invalid.");
    }

    switch (pi) {
    case 0:
      return getR0();
    case 1:
      return getA();
    case 2:
      return getB();
    case 3:
      return getC();
    case 4:
      return getD();
    default:
      throw std::runtime_error("Default switch case should not be reachable.");
    }
  }

  bool NormalizedSteinhartHartModelThermistor::isValidResistance(Resistance r) const {
    return isValidResistance(r);
  }

  /*
   * The (simplified) normalized Steinhart-Hart equation:
   *  1/T = A + B * ln(R/R0) + C * ln³(R/R0)
   */
  Temperature NormalizedSteinhartHartModelThermistor::getTemperature(Resistance r) const {
    double R0 = getR0();
    double A = getA();
    double B = getB();
    double C = getC();
    double D = getD();
  
    double lnR = log(r/R0);
    double Tinv = A + B * lnR + C * pow(lnR, 2) + D * pow(lnR, 3);
    double T = 1.0/Tinv;

    return T;
  }


  bool NormalizedSteinhartHartModelThermistor::isValidTemperature(Temperature t) const {
    return isValidKelvinTemperature(t);
  }
  Resistance NormalizedSteinhartHartModelThermistor::getResistance(Temperature t) const {
    double R0 = getR0();
    double A = getA();
    double B = getB();
    double C = getC();
    double D = getD();

    double x = C / D;
    double y = B / D;
    double z = (A - 1.0/t) / D;

    double p = y - pow(x,2)/3;
    double q = 2.0/27.0 * pow(x,3) - x*y/3 + z;

    double s = sqrt( pow(q,2)/4 + pow(p,3)/27 );
    double u = pow( -q/2 + s, 1.0/3.0 );
    double v = -pow( q/2 + s, 1.0/3.0 );

    double r = exp( u + v - x/3 );
    if (r == NAN) {
      return NAN;
    }

    return r * R0;
  }
  /*Resistance NormalizedSteinhartHartModelThermistor::getResistance(Temperature t) const {
    printf("Finding resistance for t = %.7f:\n", t);
    if (t < 273.15) {
      printf("\tNegative temperature!\n");
    }
    const double R0 = getR0();
    gsl_root_fsolver* s = gsl_root_fsolver_alloc(gsl_root_fsolver_falsepos);

    double r = 0;
    double x_lo = 0.000000001, x_hi = 1000000000000000000;
    struct root_finding_params params = {t, const_cast<NormalizedSteinhartHartModelThermistor*>(this)};

    gsl_function F;
    F.function = &getTemperatureStatic;
    F.params = &params;

    double y_lo = getTemperatureStatic(x_lo, &params);
    double y_hi = getTemperatureStatic(x_hi, &params);
    if ((y_lo < 0 && y_hi < 0) || (y_lo >= 0 && y_hi >= 0)) {
      // Same sign!
      printf("\tSame sign: x_lo=%.7f->y_lo=%.7f   x_hi=%.7f->y_hi=%.7f\n", x_lo, x_hi, y_lo, y_hi);
      return NAN;
    }
    gsl_root_fsolver_set(s, &F, x_lo, x_hi);
    

    int status;
    int iter = 0;
    const int max_iter = 100;
    do {
      iter++;
      status = gsl_root_fsolver_iterate(s);
      r = gsl_root_fsolver_root(s);
      x_lo = gsl_root_fsolver_x_lower(s);
      x_hi = gsl_root_fsolver_x_upper(s);
      status = gsl_root_test_interval(x_lo, x_hi, 0, 0.000001);
     
      if (status == GSL_SUCCESS)
	printf ("Converged:\n");
      
      printf("%5d [%.7f, %.7f] %.7f %.7f\n",
	      iter, x_lo, x_hi, r, x_hi - x_lo);
    }
    while (status == GSL_CONTINUE && iter < max_iter);

    gsl_root_fsolver_free(s);
    return r;
    }*/

  double NormalizedSteinhartHartModelThermistor::getTemperatureStatic(double r, void *data)
  {
    struct root_finding_params* params = (struct root_finding_params*)data;
    return params->thermistor->getTemperature(r) - params->t;
  }
 
  void NormalizedSteinhartHartModelThermistor::accept(ThermistorVisitor* visitor) const
  {
    visitor->visit(this);
  }

  std::ostream &NormalizedSteinhartHartModelThermistor::print(std::ostream &stream) const {
    return stream << "Thermistor with normalized Steinhart-Hart parameters:" <<
      " R0=" << boost::lexical_cast<std::string>(getR0()) <<
      " A=" << boost::lexical_cast<std::string>(getA()) <<
      " B=" << boost::lexical_cast<std::string>(getB()) <<
      " C=" << boost::lexical_cast<std::string>(getC()) <<
      " D=" << boost::lexical_cast<std::string>(getD());
  }

}

